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ABSTRACT 

We have conducted a near-infrared monitoring campaign at the UK InfraRed Telescope 
(UKIRT), of the Local Group spiral galaxy M33 (Triangulum). The main aim was to 
identify stars in the very final stage of their evolution, and for which the luminosity 
is more directly related to the birth mass than the more numerous less-evolved giant 
stars that continue to increase in luminosity. In this second paper of the series, we 
construct the birth mass function and hence derive the star formation history. The 
star formation rate has varied between ~ 0.002 and 0.007 Mq yr~^ kpc~^. We give 
evidence of two epochs of a star formation rate enhanced by a factor of a few - one 
that happened > 6 Gyr ago and produced > 80% of the total mass in stars, and one 
around 250 Myr ago that lasted ^ 200 Myr and formed < 6% of the mass in stars. 
We construct radial and azimuthal distributions in the image plane and in the galaxy 
plane for populations associated with old first-ascent red giant branch (RGB) stars, 
intermediate-age Asymptotic Giant Branch (AGB) stars and young (massive) blue and 
red supergiants. We find that the RGB stars follow a spheroidal distribution, while 
younger stars follow a flat-disc distribution. The intermediate-age population displays 
signs of a pseudo-bulge or possibly a bar. The inner spiral arm pattern as recorded 
in mid-19'^-century drawings is confirmed. We interpret our findings as evidence for 
an old, pressure-supported component and a younger disc formed 6 Gyr ago, with an 
accretion event occurring 250 Myr ago giving rise to the compact nucleus in M33. Our 
study provides support for recent Padova stellar evolution models except that super- 
AGB stars likely reach low temperatures and thus high mass-loss rates, supporting the 
super- AGB nature of the progenitors of dust-enshrouded supernovae such as SN 2008S. 

Key words: stars: evolution - stars: luminosity function, mass function - galax- 
ies: individual: MSB - galaxies: star formation ~ galaxies: stellar content - galaxies: 
structure 



1 INTRODUCTION 

Spiral galaxies are among the more massive galajcies in the 
Universe, that in the A-Cold Dark Matter paradigm are be- 
lieved to have formed from the agglomeration of smaller 
building blocks. They have gaseous, rotating discs, where 
spiral arms and rotational shear set the conditions within 
which star formation continues at a moderate rate in present 
times. In the central regions of spiral galaxies, more ellip- 
soidal, largely stellar and pressure-supported components 
exist - nuclear star clusters, bulges, and "bars" , but the rel- 
ative importance of these components varies among galaxies 
of this Hubble class. What determines the presence or ab- 
sence of these structures, and indeed the spiral arm pattern 
itself, is not understood. 

The Local Group galaxy Triangulum (Hodierna 1654) 



- hereafter referred to as M 33 (Messier 1771) - offers us a 
unique opportunity to study a spiral galaxy up close, and 
in particular to learn more about the structure and evolu- 
tion of the central regions of such galaxy, that in our own 
Milky Way are heavily obscured by the intervening dusty 
Disc (van Loon et al. 2003; Benjamin et al. 2005). Our view- 
ing angle with respect to the M 33 disc is more favourable 
(56-57° - Zaritsky, Elston & Hill 1989; Deul & van der Hulst 
1987) than that of the larger M31 (Andromeda), whilst the 
distance to M 33 is not much different from that to M 31 
(/i — 24.9 mag - Bonanos et al. 2006). The stellar content 
of M 33 was summarised beautifully in the review by van 
den Bergh (1991). The present-day star formation rate in 
the central regions of M 33 has been measured on the basis 
of optical photometry of massive main-sequence stars and 
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Ha luminosity to be ~ 0.005-0.1 Mq yr"^ kpc^^ (Wilson, 
Scoville & Rice 1991). Minniti, Olszewski & Rieke (1993) de- 
duced from their infrared (IR) observations (in the H-band 
at 1.6 pm) that a distinct star formation episode must have 
occurred in the "bulge" of M33 within the last Gyr, unre- 
lated to activity within the disc - they also showed that the 
bulge comprises also intermediate-age stars in contrast to 
the old, metal-poor halo of M33. The structure of the cen- 
tral region was investigated further at IR wavelengths by 
Regan & Vogel (1994) and Stephens & Frogel (2002); the 
central nucleus had been described in detail also by Lauer 
et al. (1998). 

Galactic evolution is driven at the end-points of stel- 
lar evolution, where copious mass loss returns chemically- 
enriched and sometimes dusty matter back to the interstel- 
lar medium (ISM); the stellar winds of evolved stars and the 
violent deaths of the most massive stars also inject energy 
and momentum into the ISM, generating turbulence and 
galactic fountains when superbubbles pop as they reach the 
"surface" of the galactic disc. The evolved stars are also ex- 
cellent tracers, not just of the feedback processes, but also of 
the underlying populations, that were formed from millions 
to billions of years prior to their appearance. The evolved 
phases of evolution generally represent the most luminous, 
and often the coolest, making evolved stars brilliant beacons 
at IR wavelengths, where it is also easier to see them deep 
inside galaxies as dust is more tranpsarent at those longer 
wavelengths than in the optical and ultraviolet where their 
main-sequence progenitors shine. The final stages of stellar 
evolution of stars with main-sequence masses up to M ~ 30 
M0 - Asymptotic Giant Branch (AGB) stars and red su- 
pergiants - are characterised by strong radial pulsations 
of the cool atmospheric layers, rendering them identifiable 
as long-period variables (LPVs) in photometric monitoring 
campaigns spanning months to years (e.g., Whitelock, Feast 
& Gatchpole 1991; Wood 2000; Ita et al. 2004a,b). 

The main objectives of the project are: to construct the 
mass function of LPVs and derive from this the star forma- 
tion history in M 33; to correlate spatial distributions of the 
LPVs of different mass with galactic structures (spheroid, 
disc and spiral arm components); to measure the rate at 
which dust is produced and fed into the ISM; to estab- 
lish correlations between the dust production rate, luminos- 
ity, and amplitude of an LPV; and to compare the in situ 
dust replenishment with the amount of pre-existing dust 
(see Javadi, van Loon & Mirtorabi 2011b). This is Paper 
II in the series, describing the galactic structure and star 
formation history in the inner square kpc. Paper I in the 
series presented the photometric catalogue of stars in the 
inner square kpc (Javadi, van Loon and Mirtorabi 2011a). 
Subsequent papers in the series will discuss the mass-loss 
mechanism and dust production rate (Paper III), and the 
extension to a nearly square degree area covering much of 
the M33 optical disc (Paper IV). 



2 THE INPUT DATA AND MODELS 

In this section we briefly summarise the salient features of 
the catalogue of variable stars that we produced (Paper I), 
and of the stellar evolution models that we use here. 




(J-K) 



Figure 1. Colour— magnitude diagram of our UKIRT catalogue; 
variable stars are plotted in green. Representative photometric 
errorbars are plotted. Overlain are isochrones (Marigo et al. 2008), 
labelled by their logarithmic ages - these include circumstellar 
reddening, showing up as large excursions towards red colours; 
carbonaceous dust (e.g., at logt = 9) and oxygeneous dust (e.g., 
at logt = 10) have different reddening slopes. The vertical line at 
J—K= 1.5 mag indicates the colour criterion redwards of which a 
correction is applied for this reddening. The thick black-and-white 
lines demarcate the bulk of the populations of massive, AGB and 
RGB stars (separated by small buffers in K-band magnitude to 
avoid cross-contamination - see Section 4.2). 



2.1 The catalogue of variable stars that we use 

In Paper I we described the search for large- amplitude, long- 
period variable stars in the central 4' x 4' (square kpc) of 
M 33, with the UIST imager on the United Kingdom IR Tele- 
scope (UKIRT) on Mauna Kea, Hawai'i. The monitoring was 
done predominantly in the K-band (around a wavelength of 
A = 2.20 /xm), with additional observations in the J- and H- 
bands (A = 1.28 and 1.65 ^im, respectively) mostly for the 
purpose of obtaining colour information. Observations were 
carried out if the near-IR stellar images were < 0.8". Each 
quadrant of the field was observed at least 11 times spread 
over the period 2003-2007, with stars in overlap regions be- 
ing observed more often. 

Photometry was obtained via fitting of the empirically- 
derived Point Spread Function to all stars detected above a 
threshold in terms of the local sky background noise. Hence 
a catalogue of 18,398 stars was produced. Fig. 1 shows these 
stars in a colour-magnitude diagram. The catalogue reaches 
a completeness level of > 90% above the RGB tip around 
7^ ~ 18 mag. For this bright magnitude regime, the photo- 
metric accuracy is < 10%, and blending is a concern in only 
~ 1% of cases. Contamination of the sample by stars in the 
Galactic foreground was demonstrated to be negligible. 

The variability search was performed using the new- 
trial routine (Stetson 1993) on the basis of the "L" vari- 
ability index (Stetson 1996), which includes an assessment 
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of the variations in the measurements in comparison to the 
measurement uncertainties, both in terms of their ratio and 
distribution. Hence 812 variables were found, with estimated 
K-band ampUtudes in the range Ay, ~ 0.2-2 mag (typi- 
cally Ak ^ 0.6 mag). Their distribution over magnitude 
and colour is shown in Fig. 1 (in green). 

2.2 The stellar evolution models that we use 

The best available theoretical models for our purpose are 
currently provided by the Padova group (Marigo et al. 2008), 
for the following reasons: 

• They span an adequate range in birth mass (0.8 < M < 

30 M©), combining their own models for intermediate-mass 
stars {M < 7 Mq) with Padova models for more massive 
stars (M > 7 Mq; BcrtcUi et al. 1994); 

• They are computed all the way through the thermal 
pulsing AGB until they enter the post-AGB phase, in a 
manner that is consistent with the computation of the pre- 
ceding evolutionary stages. Crucially, it includes the third 
dredge-up mixing of the stellar mantle as a result of the 
helium-burning pulses, as well as the enhanced luminosity of 
massive AGB stars undergoing Hot Bottom Burning (HBB 
- Iben & Renzini 1983); 

• They include molecular opacities that are important for 
the cool atmospheres of red giants, allowing to describe the 
transformation from oxygen-dominated M-type AGB stars 
to carbon stars in the approximate birth-mass range M ^ 
1.5-4 Mq (cf. Girardi & Marigo 2007); 

• They include predictions for dust production in the 
winds of LPVs and the associated reddening; 

• They include predictions for the radial pulsation; 

• They have been carefully transformed to various optical 
and IR photometric systems; 

• They arc available via a user-friendly internet interface. 

We encourage other groups to apply the same principles to 
their models, so that their fundamental differences can be 
explored empirically to the greatest effect. While we refer to 
Marigo et al. (2008) and Marigo & Girardi (2007) for the 
technical details pertaining to their models, and we have no 
way of replacing them by different models, we point out a 
few essential aspects of this set of Padova models below. 

The Padova models do not include the effects of rota- 
tion or magnetic fields on mixing within the stellar interior. 
Especially in massive stars these can have important con- 
sequences for their luminosities and lifetimes, making them 
brighter but live longer and evolve to cooler red supergiant 
stages (cf. Meynet & Maeder 2000; Heger, Woosley & Spruit 
2005). 

The Padova models employ convectivc overshoot (sec 
Marigo & Girardi 2007), which prolongs the lifetime of a 
star. Hence stars of a given observed luminosity will be 
linked to older ages and lower birth masses than in other, 
classical models that do not allow for convectivc overshoot. 
This is mainly an issue for intermediate-age populations. 

The thermal pulsing AGB is not computed by solv- 
ing the equations governing the stellar internal structure, 
but in a synthetic way by adopting paranieterisations which 
include the efficiency and threshold temperature at which 
third dredge-up sets in. These parameters have been cali- 
brated empirically against observations, for instance of the 



carbon star luminosity function and star clusters in the Mag- 
ellanic Clouds (Marigo & Girardi 2007). Consequently, it is 
possible that these calibrations are not a perfect match to 
populations in galaxies other than the Magellanic Clouds. 
While stars in M 33 are not vastly different from those in 
the LMC, stars in the nuclear region of M33 are probably 
more similar to solar-metallicity stars in the Milky Way. 

More critically important is to note that the published 
isochrones do not include the excursions in luminosity and 
temperature resulting from the thermal pulses, but only 
record the evolution through the hydrogen-shell burning 
phases. While the latter dominates also on the TP-AGB, 
stars undergoing the thermal-pulse excursions may erro- 
neously be assigned to lower or higher birth masses. This 
ultimately reduces the time resolution and contrast of the 
derived star formation history somewhat. 

The radial pulsation properties (period and mode) are 
computed as a function of luminosity. Stars arc expected to 
start pulsating in the first overtone as soon as the thermal 
pulses start. The luminosity at which the transition to the 
fundamental mode occurs is based upon the linear pulsation 
models from Ostlie & Cox (1986), while periods are com- 
puted from the mass and radius based upon models from 
Fox & Wood (1982) - with some modifications to detail as 
described in Marigo & Girardi (2007; their Section 2.7) and 
Marigo et al. (2008; their Appendix). For the purpose of our 
work the main important feature is the lifetime of the pul- 
sation phase, not so rrmch the mode or period. The revisions 
made by Marigo et al. (2008) reduced the pulsation lifetime 
of the most massive AGB stars by a factor of 2-5, which in 
our analysis would lead to an increase in the associated star 
formation rate by the same amount (see Section 3.2). 

Spectral templates used to describe the photospheric 
emission of the coolest giants come from empirical libraries 
for M-type stars (Fluks et al. 1994) or synthetic spectra for 
carbon stars (Loidl, Lancon & j0rgcnsen 2001). These may 
not always be entirely appropriate, for instance at a different 
metallicity, carbon-to-oxygen ratio or dynamical state of the 
pulsating atmosphere. But the uncertainty in the bolomctric 
corrections is modest in the near-IR, likely < 0.1 mag at K, 
and thus not a serious concern for our project. 

Mass loss truncates the evolution of AGB stars, deplet- 
ing their mantles and thereby avoiding the explosion of AGB 
stars of M > 1.4 M©. This becomes important when the 
mass-loss rates exceed the nuclear burning rate (van Loon 
et al. 1999). Mass loss also reduces the effective temperature 
and increases the pulsation period, as the less-massive man- 
tle inflates. Because all this happens quite suddenly, it is as 
important to accurately predict when this happens, as it is 
to accurately predict the rate of mass loss that is attained. 
The formalism adopted in the Padova models is based upon 
parameterisation with respect to pulsation period (cf. Vas- 
sialidis & Wood 1993). While this has observational sup- 
port it does mean that the accuracy of the predicted mass 
loss relies critically on the model's accuracy with which it 
predicts the pulsation properties. Alternative prescriptions 
for mass loss exist that are defined in terms of luminosity 
and temperature (van Loon et al. 2005). Mass loss is not as 
dramatic an issue for massive stars, which do explode be- 
fore shedding their mantle through stellar mass loss, though 
mass loss could have an effect on the time these stars spend 
as a red supergiant. 
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Mass loss from cool stars also afltects the photometric 
appearance of the star, as dust forming in the wind gives 
rise to selective extinction and thus diminishing brightness 
and reddening at optical/near-IR wavelengths. The effect at 
near-IR wavelengths only becomes noticeable in the case of 
very dense circumstellar envelopes. The Padova models in- 
corporate these effects in a two-step procedure: (1) predict 
the radial density of the dust envelope, and (2) compute 
the resulting extinction of starlight at short wavelengths, 
and re-emission by the grains at long wavelengths. The first 
step is based mainly on the predictions from Ferrarotti & 
Gail (2006) whilst the second step follows standard radia^ 
five transfer computations as in Grocncwcgcn (2006). The 
adopted wind structure and dust content have observational 
support (e.g., Marshall et al. 2004) but there can be many 
reasons for deviations from this standard, simplistic pic- 
ture. Likewise, the optical properties of the grains that are 
adopted in the radiative transfer computations are reason- 
able (typical astronomical silicates or aluminium oxides for 
M-type stars and a mixture of amorphous carbon and sili- 
cate carbide for carbon stars) but not without uncertainties. 
For the purpose of the analysis presented here, the most 
important uncertainty is in the direction of the reddening 
vector in the near-IR colour-magnitude diagram. This is 
explored in some detail in Sections 3.1.1 and 4.1.1. 

An important remaining uncertainty in the Padova 
models concerns the evolution of super-AGB stars, with 
birth-masses M ~ 5-10 Mq (Siess 2007), which in their 
models are not computed all the way through the thermal 
pulsing phase and thus appear to terminate their evolution 
prematurely. We explore the consequences in this paper and 
offer some guidance to the models. 

Fig. 1 includes several representative isochroncs, which 
show good correspondence to the observed sequences of 
stars. The main regions in this diagram where massive stars, 
AGB stars and RGB stars arc found are demarcated (see 
Section 4.2). Red AGB stars or massive stars, at J — K > 1.5 
mag, are almost exclusively LPVs affected by circumstellar 
reddening; interstellar reddening was proven to be unim- 
portant (see Paper I). The reddened stars form a minority 
among the sample of LPVs, but it is encouraging to sec that 
the isochrones cover the regime of reddened LPVs fairly well. 



3 FROM THE BRIGHTNESSES OF VARIABLE 
STARS TO A STAR FORMATION HISTORY 

The star formation history (SFH) is described by the star 
formation rate (SFR), ^, as a function of lapsed time, t. 
This function quantifies how many solar masses of gas are 

converted into stars per year, and we here express it per 
unit area projected on the sky. The combined mass of stars 
created between times t and t + dt is: 



dM{t) = ^ (t) dt. 

This corresponds to a number, N, of stars formed: 
/.!lir/iMF(m)dm 



dN{t) = dM{t)- 



f /imp (m)m dm ' 
where /imp is the initial mass function, defined by: 
/imp = Am'", 



(1) 



(2) 



(3) 



where A is the normalization constant and a depends on the 
mass range, following Kroupa (2001): 



{+0.3 ± 0.7 for min < m/Mo < 0.08 
-1-1.3 ± 0.5 for 0.08 < m/Mo < 0.50 
-1-2.3 ± 0.3 for 0.50 < m/M® < max 



(4) 



The minimum and maximum of the stellar mass range are 
adopted to be 0.02 and 200 M®, respectively. Reasonable 
changes in these values will result in changes in the star 
formation rate by hardly more than a factor two. 

The question is how many of these stars, n, are variable 
stars around time t which is when we observe them. Suppose 
stars with mass between m{t) and m(t + dt) satisfy this 
condition, then the number of variable stars created between 
times t and t + dt is: 



dn(i) = dN{t) 



m(t) 



/iMp(m) dm 



/rir/iMp(m)dm 
Combining the above equations we have: 



dn{t) = £.{t)dt 



Jm(t) /iMp(m)dm 
f /imp (m)m dm 

Jmin ^ ' 



(5) 



(6) 



When wc determine ^{t) over an age bin dt, the number 
of variable stars observed in that age bin, dn' , depends on 
the duration of variability, 6t, so: 



dn'(t) = dn(t)^. 



(7) 



Inverting the above equations we obtain a relation between 
the observed number of variable stars in a certain age bin, 
and the star formation rate that long ago: 



m = 



, , ^ f max « / \ 1 

(t) Jj„i„ fiMF{rn)rndm 



St 



I 

J r 



m(t+dt) 
>(t) 



/imp (m) dm 



(8) 



Note that the value of the normalization constant A in Eq. 
(3) does not matter here. 



3.1 The relation between birtli mass and K-band 
magnitude of large-amplitude VEtriable stars 

LPVs, which dominate the large-amplitude variables iden- 
tified in our IR monitoring programme, have reached the 
very final stages of their evolution, and their brightness can 
thus be translated into their mass at birth by employing 
theoretical evolutionary tracks, or in this case just as well 
by isochrones. Inspection of the isochrones confirms that 
stars are at their brightest in the K-band (A ~ 2 fim) when 
they also develop the pulsations that give rise to their large- 
amplitude variability. 

Hence we constructed the mass-luminosity relation for 
the K-band as depicted in Fig. 2 (low- to intermediate-mass 
range) and Fig. 3 (which includes the most massive stars), 
for four different values of the overall metallicity, from super- 
solar - as in massive elliptical galaxies and sub-populations 
in the bulges of massive spiral galaxies such as the Milky 
Way - to sub-solar values appropriate for the Largo and 
Small Magellanic Clouds. There are small but noticeable 
differences. The central region of M33 is characterised by 
approximately solar-metallicity ISM (Rosolowsky & Simon 
2008; Magrini et al. 2009) and young and intermediate-age 
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Figure 2. Mass-Luminosity relation for Z = 0.019, 0.015, 0.008 
and 0.004, for the K-band magnitude range corresponding to the 
low- to intermediate-mass range. Sohd lines are the best linear 
spline fits, while the dotted lines are the fit for Z = 0.015 - which 
we apply for the variables stars of M33 galaxy - for comparison. 



Figure 3. Mass-Luminosity relation for Z = 0.019, 0.015, 0.008 
and 0.004. This graph shows the relation for the full range of 
K-band magnitudes (dots). The solid lines are the linear spline 
fits, for the case in which the function is interpolated across 
the super-AGB phase to massive red supergiants, i.e. for 0.7 < 
log(M/MQ) < 1.2-1.3. 



stellar populations (Massey et al. 1996; Minniti et al. 1993), 
but also given that older stars will tend to have lower metal- 
licities we err on the side of caution and adopt the Z = 0.015 
result, which is either consistent with solar metallicity or 
slightly sub-solar metallicity depending on the exact metal- 
licity adopted for the Sun {Z — 0.015 is close to the value 
for the Sun determined in recent years - cf. Asplund et al. 
2009). 

The mass-luminosity relations show small excursions 
in places, e.g., for 16 < K < 16.4 mag sd, Z — 0.019, 
which could be related to differences in the bolometric cor- 
rections of M-type and carbon stars. However, some inaccu- 
racies result from limited accuracy in the models as well as 
the time-sampling of the evolutionary tracks and isochrones, 
and where there is doubt we interpolate over these small ex- 
cursions by means of a spline fit, rendering a well-behaved 
function. 

For 0.7 < log(M/M0) < 1.2-1.3, the isochrones suggest 
an excursion towards fainter K-band magnitudes (Fig. 3), 
but as noted earlier the models were not computed through 
to the end of the evolution of these super-AGB stars (and 
possibly stars at the lower end of the mass range for red 
supergiants). In our analysis, we follow two different ap- 
proaches, and investigate the impact on the derived star 
formation history: 

1. We interpolate the Mass-Luminosity realtion (and the 
Mass-Pulsation duration relation, see further) over this 
mass range to connect the final luminosities of AGB stars 
to those of massive red supergiants in a monotonic manner. 
Stars with log(M/MQ) > 1.5 do not appear to become red 
supergiants and they are therefore not expected to exhibit 
periodic large-amplitude variability. This is our preferred so- 
lution, which we vindicate in Section 4.1.5; 



2. We accept the jump to fainter magnitudes from AGB to 
super-AGB stars. This is an extreme measure, as we believe 
super-AGB stars should reach high luminosities as they un- 
dergo HBB, though more massive red supergiants do not. 
We postpone following this choice to Section 4.1.5, where 
we show that it leads to physically unrealistic results. 

The coefficients of the best fitting function are obtained by 
minimisation of with the Iraf task nfitId except where 
it was deemed more appropriate to be guided by two end- 
points of isochrones, in which case no estimates of uncer- 
tainty are attached. For case (1) the function is described in 
Table 1. 

3.1.1 Correction for extinction by dust 

By determining the Mass-Luminosity relation we now can 
convert the K-band magnitudes of the variable stars to their 
masses, and derive the present-day mass function of variable 
stars. However, the K-band magnitude of a star seen through 
a high dust column density is diminished, and should be 
corrected to the value on the Mass-Luminosity relation. Al- 
though the isochrones of Marigo et al. (2008) include the 
effects of circumstellar dust formation, these are very uncer- 
tain and we prefer a simple de-reddening procedure. This 
would also naturally account for any additional interstellar 
extinction, at least to some extenlp. The de-reddening cor- 

^ Interstellar dust will have somewhat different optical properties 
from those of circumstellar dust. Then again, interstellar extinc- 
tion is not generally observed to be a major factor in this region 
of M33, whereas circumstellar dust is evident in a considerable 
fraction of the highly-evolved AGB stars - see Paper I. 
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Table 1. Mass— Luminosity relation, logM [Mg] = aK + b, where 
K is the K-band magnitude on the 2MASS system, for four values 
of metallicity, for the case in which the function is interpolated 
across the super- AGB phase to massive red supergiants. 



Z = 0.019 


a 


b 


validity range 


-1.627 

-0.458 
-0.188 ±0.053 
-0.094 ± 0.045 
-0.532 ± 0.033 
-0.225 ± 0.029 
-0.143 ± 0.041 


23.3812 
6.5002 
3.492 ± 0.081 
2.027 ±0.071 
9.049 ± 0.053 
3.951 ±0.045 
2.541 ±0.061 


K < 13.590 
13.590 < K < 14.875 
14.875 < K < 15.462 
15.462 < K < 16.049 
16.049 <K < 16.635 
16.635 <K < 17.222 

K > 17.222 


Z = 0.015 


a 


b 


validity range 


-2.198 

-0.535 
-0.213 ±0.027 
-0.088 ±0.018 
-0.222 ± 0.021 
-0.582 ± 0.012 
-0.147 ±0.013 
-0.173 ± 0.022 


31.3051 
8.6253 
3.844 ±0.040 
1.924 ±0.023 
4.056 ± 0.024 
9.937 ±0.014 
2.620 ±0.016 
3.060 ±0.032 


K < 13.633 
13.633 < K < 14.854 
14.854 <K < 15.347 
15.347 < K < 15.840 
15.840 < K < 16.332 
16.332 <K < 16.825 
16.825 < K < 17.318 

K > 17.318 


Z = 0.008 


a 


b 


validity range 


-0.840 

-0.589 
-0.188 ±0.028 
-0.142 ±0.027 
-0.188 ±0.022 
-0.501 ± 0.016 
-0.248 ± 0.018 
-0.128 ± 0.025 


12.840 
9.391 
3.465 ± 0.043 
2.758 ±0.041 
3.487 ±0.036 
8.183 ±0.024 
4.335 ± 0.026 
2.257 ±0.035 


K < 13.732 
13.732 < K < 14.787 
14.787 < K < 15.300 
15.300 < K < 15.813 
15.813 < K < 16.327 
16.327 <K < 16.840 
16.840 < K < 17.353 

K > 17.353 


Z = 0.004 


a 


b 


validity range 



-0, 
-0.209 
-0.240 
-0.050 
-0.131 
-0.243 
-0.714 
-0.109 
-0.153 



708 

± 0.051 
± 0.054 
± 0.060 
± 0.054 
± 0.041 
± 0.039 
± 0.043 
± 0.046 



11.1268 
3.783 ± 0.077 
4.244 ± 0.078 
0.565 ±0.091 
2.583 ±0.083 
4.409 ± 0.057 
12.310 ±0.055 
1.932 ±0.070 
2.690 ±0.067 



K 

14.734 
15.140 
15.545 
15.951 
16.356 
16.762 
17.167 
K 



< 14.734 

< K < 15.140 

< K < 15.545 
<K < 15.951 

< K < 16.356 

< K < 16.762 
<K < 17.167 

< K < 17.573 
> 17.573 



correction to those variable stars that have J — K > 1.5 
mag to allow some deviations from the unreddened locus as 
a result of photometric inaccuracies and variability. Hence 
the correction equation is: 



Ka = K + a{1.2b - {J - K)) 



(9) 



If only K- and H-band data are available but no J-band data 
then we apply the following correction derived from the K 
versus H-K diagram: 



Ko = K + a{G.2,-{H~K)) 



(10) 



In the case where only K-band but no J- nor H-band data 
are available, we adopt J—K = 4 mag (roughyl equivalent to 
where we start losing stars too red for our survey depth), and 
use Eq. (9). The constants in these equations are a = 1.06 
for carbon stars and a — l.Q for oxygen-rich star^. 

It is thus important to know which stars have carbona- 
ceous dust and which ones have oxygeneous dust. As de- 
scribed in Paper I, the only catalogue available that uses an 
identification technique based on specific spectral features 
is the one published by Rowe et al. (2005). Then, follow- 
ing the recommendation by those auhtors, we set limits on 
the photometric errors in their catalogue of < 0.05 mag, 
resulting in just three variable carbon stars. This is insuf- 
ficient to decide which red stars in our catalogue are likely 
to be carbon stars. Setting no limit on errors resulted in 
five variable carbon stars, of which three have J — K > 1.5 
mag. However, we expect from theory and observations in 
the LMC (Groenewegen & de Jong 1993; van Loon, Mar- 
shall & Zijlstra 2005; Girardi & Marigo 2007) that solar and 
slightly sub-solar metallicity carbon stars arise from stars 
in a birth mass range of ~ 1.5-4 Mq: 3'^'* dredge-up is not 
strong enough to turn lower-mass AGB stars into carbon 
stars, while HBB prevents carbon to enrich the surface of 
more massive AGB stars. Hence we first apply the carbon- 
dust correction to all reddened stars - if the derived mass is 
in the 1.5-4 M© range then we accept that star as a carbon 
star, else we apply the oxygeneous dust correction. 

In Fig. 4 we show the resulting PMF for the variable 
stars in the central square kpc of M 33. As indicated in Paper 
I, nine of the large-amplitude variable stars show indications 
of blending so we removed those stars from our anlysis and 
we derive the star formation history using the remaining 
803 variables. We applied the de-reddening correction as de- 
scribed above, as well as only a carbonaceous dust correction 
and only an oxygeneous dust correction: the differences are 
slight and mostly reflected in an sin 50% scatter in the ~ 2 
M0 region. 



rection is related to the colour and magnitude of a star. In 
Fig. 1 we show the Marigo et al. (2008) isochrones for a 
range of ages. It immediately becomes clear that the red- 
dening for carbon stars (e.g., at t = 1 Gyr i.e. \ogt — 9) 
differs from that of oxygen-rich stars (e.g., at t = 10 Gyr i.e. 
logt = 10), the former more quickly reddening compared to 
the extinction in the K-band. The average slope of the red- 
dening of carbon stars is 0.52 mag of K-band extinction per 
magnitude oi J — K reddening, whilst it is 0.72 mag mag~^ 
for oxygen-rich stars. The locus to which we translate the 
reddened star is at J — K = 1.25 mag, but we apply the 



3.2 The mass-dependent duration of 
large-amplitude variability 

Eq. (1) contains a correction factor for the duration in which 
a star is seen to pulsate with large amplitude. For instance, 
massive stars evolve faster and are captured in their pulsat- 
ing phase only briefly, whereas lower-mass stars spend longer 
in that phase and are thus more likely and more numerous 



^ In the remainder we assume all photometry has been corrected 
for extinction, and we do not explicitly use the subscript "0" . 
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Figure 4. Present-day mass function of large-amplitude variable 
stars in the central square kpc of M33. The solid black line is 
the histogram obtained by applying the de-reddening correction 
dependent on the mass of the star, the dotted red and dashed 
blue lines are the results obtained when applying only the car- 
bonaceous dust correction or only the oxygeneous dust correction, 
respectively. 



Table 2. Mass-Pulsation relation, log(<5t/t) = D + 

X^Ll "^^P [^©] ~ *i)^/(2c^)) , where St is the pulsa- 

tion duration and t the age of the star, for four values of metal- 
licity, for the case in which the function is interpolated across the 
super-AGB phase to massive red supergiants. 



Z = 0.019 


D 


i 


a 


6 


c 


-4.42 


1 


3.04 


1.271 


0.277 




2 


0.74 


0.510 


0.068 




3 


1.35 


0.478 


0.331 


Z = 0.015 


D 


i 


a 


b 


c 


-3.63 


1 


2.18 


1.238 


0.238 




2 


1.36 


0.514 


0.127 




3 


0.30 


0.228 


0.084 


Z = 0.008 


D 


i 


a 


b 


c 


-3.96 


1 


2.34 


1.281 


0.378 




2 


1.32 


0.460 


0.165 




3 


0.38 


0.145 


0.067 


Z = 0.004 


D 


i 


a 


b 


c 


-4.00 


1 


2.29 


1.217 


0.408 




2 


0.84 


0.524 


0.093 




3 


0.87 


0.206 


0.088 



to appear in variability surveys - aside from the eflFects of 
the IMF, and of the SFH that we are after. 

To derive this correction we followed the predictions 
from the Marigo et al. (2008) isochrones, so as to be consis- 
tent within our analysis. The isochrone tables indicate the 
phase of strong radial pulsation, and we take this entire du- 
ration as the time during which the star would have been 
identified in our survey (Paper I) as a variable star. We de- 
rived a function of mass for this duration of pulsation by fit- 
ting a multiple- Gaussian function to these isochrone values. 
As we investigate two cases - one in which we interpolate 
between the most massive AGB stars and the more mas- 
sive red supergiants, and another in which we accept the 
regression in K-band brightness over the super-AGB star 
mass regime - we also take these two approaches with re- 
gard to the pulsation-duration versus mass function (the 
Mass-Pulsation relation): if the super-AGB stars do evolve 
to higher luminosities and cooler temperatures then they 
will also be more likely to, and longer, pulsate. 

Hence, for case (1) we fitted three Gaussian funtions to 
the Mass-Pulsation diagram and interpolated over the same 
mass range as in the Mass-Luminosity diagram. The result 
is tabulated in Table 2 and shown in Fig. 5 (top panel). For 
case (2), we fit two independent sets of Gaussian functions 
- one for the low-mass stars and one for the high-mass stars 
(see Section 4.1.5). 



-1 t 




• z = 


0.019 


■ z = 


0.015 


z = 


0.008 


« z = 


0.004 







• z = 


0,019 






■ Z = 


0.015 






Z = 


0.008 


2 




«Z = 


0.004 




log(M/Mj 

Figure 5. Mass-Pulsation relation. The points show the ratio of 
pulsation duration to age as tabulated in the Marigo et al. (2008) 
isochrones; the solid lines are our multiple-Gaussian fits to these 
data: interpolating over the super-AGB regime (top panel, for 
four choices of metallicity Z) and accepting the jump in K-band 
brightness between the AGB and super-AGB (bottom panel). 
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Table 3. Mass-Age relation, \ogt [yr] = a logM [Mg] 



0.5 1 

log(M/Mj 

Figure 6. (Birth) Mass-Age relation for AGB stars and red su- 
pergiants derived from the Marigo et al. (2008) isochrones for four 
choices of metallicity (dots, Z), along with linear spline fits. 



3.3 From a mass distribution to a star formation 
rate as a function of time 

Here we have all the relations required to calculate the SFR, 
except that we still need to relate the birth mass of a variable 
star we see now, to the time lapsed since its formation. This 
transformation is easily derived from the Marigo et al. (2008) 
isochrones. The Mass-Age relation is displayed in Fig. 6, in 
which the data points are derived from the isochrones and 
the solid line is a function we fit to the data (the coefficients 
of which are listed in Table 3). 

Now we have got all the tools ready to derive the SFIf. 
For each star, first we examine whether the star is reddened 
or not and if so, we add the de-reddening correction using the 
mass-dependent coefficients. Then, by using the coefficients 
from Table 1 for Z — 0.015, we convert the K-band mag- 
nitude to the birth mass and by using Table 2 we estimate 
the pulsation duration, the inverse of which is the weight 
assigned to the star. We do this by interpolation across the 
super- AGB regime - later (in Section 4.1.5) we investigate 
the difference obtained by following the isochrones exactly. 
Then, we convert the birth mass to lapsed time (age). The 
only remaining correction to be applied is the IMF. 

A not unimportant detail related to the presentation of 
a SFH is the way it is binned in age. The younger, massive 
variable stars are many times fewer than the old, low-mass 
variable stars and inadequate binning can either lead to spu- 
rious peaks in the SFR or mask any such real bursts. From 
a statistical point of view there is advantage is assuring that 
each bin contains the same number of stars, so the SFR val- 
ues have uniform uncertainties. To accomplish this, we first 
ordered the stars by mass, then start counting until we reach 
a given number, at which moment we start counting stars 
for the subsequent bin. Eventually, the IMF correction is 



: 0.019 



validity range 



-3.480 ± 0.008 10.082 ± 0.019 log M < 0.200 

-2.433 ±0.007 9.873 ±0.010 0.200 < log < 0.415 

-2.866 ±0.007 10.054 ±0.010 0.415 < log M < 0.631 

-2.381 ±0.007 9.750 ±0.010 0.631 < log M < 0.846 

-2.027 ±0.008 9.450 ±0.012 0.846 < log M < 1.061 

-1.463 ±0.009 8.852 ±0.013 1.061 < log M < 1.277 

-1.168 ±0.011 8.474 ±0.016 1.277 < log M < 1.492 

-0.621 ±0.014 7.600 ±0.019 logAf> 1.492 



: 0.015 



validity range 



-3.062 ±0.007 10.096 ±0.011 logM< 0.237 

-2.425 ± 0.007 9.945 ± 0.010 0.237 < log M < 0.455 

-2.824±0.007 10.127±0.010 0.455 < log M < 0.674 

-2.316 ±0.007 9.786 ±0.011 0.674 < log M < 0.892 

-1.926 ±0.008 9.438 ±0.012 0.892 < log iW < 1.110 

-1.399 ±0.009 8.853 ±0.013 1.110 < logM < 1.328 

-1.180 ±0.063 8.562 ±0.067 1.328 < logM < 1.546 

-0.625 ±0.013 7.704 ±0.018 logM > 1.546 



: 0.008 



validity range 



-3.461 ±0.008 9.976 ±0.012 logM < 0.179 

-2.347 ±0.007 9.776 ±0.011 0.179 < logM < 0.404 

-2.727 ±0.008 9.930 ±0.011 0.404 < log M < 0.628 

-2.154 ±0.008 9.570 ±0.014 0.628 < log M < 0.853 

-1.848 ±0.009 9.309 ±0.015 0.853 < log M < 1.077 

-1.398 ±0.010 8.825 ±0.015 1.077 < logM < 1.302 

-1.134 ±0.012 8.451 ±0.017 1.302 < logM < 1.526 

-0.681 ±0.015 7.790 ±0.021 logM > 1.526 



: 0.004 



validity range 



-3.416 ±0.007 9.910 ±0.011 logM < 0.161 

-2.436 ±0.007 9.752 ±0.010 0.161 < logM < 0.389 

-2.588 ±0.007 9.812 ± 0.011 0.389 < log M < 0.617 

-2.043 ± 0.007 9.475 ± 0.011 0.617 < log M < 0.845 

-1.812 ±0.008 9.280 ±0.012 0.845 < log M < 1.073 

-1.395 ±0.010 8.832 ±0.013 1.073 < logM < 1.301 

-1.036 ±0.011 8.365 ±0.016 1.301 < logM < 1.529 

-0.829 ±0.014 8.048 ±0.020 logM > 1.529 



applied based on the minimum and maximum mass of each 
bin. 



4 THE STAR FORMATION HISTORY IN M 33 

The obtained SFH for the central square kpc of M 33 is 
shown in Fig. 7. The horizontal "error bars" correspond to 
the start and end time for which the SFR was calculated in 
that bin. The oldest bin is only shown in part as it extends 
to (unrealistically) large ages. The only reason it is plotted 
is to show that, as expected, the star formation rate we 
estimate for ages exceeding the Hubble time is negligible. 
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Figure 7. The star formation liistory in tlie central square icpc 
of M 33 derived from near-IR AGB-star and red-supergiant vari- 
ables. 



We prefer not to rescale the diagrams so as not to degrade 
their visibility in the ranges where it matters. 

The SFH is characterised by a major epoch of formation 
« 4-8 Gyr ago (logt = 9.6-9.9), peaking around 6 Gyr ago 
(logt = 9.8) at a level about three times as high as during 
the subsequent couple of Gyr. This corroborates the result 
obtained by Cioni et al. (2008) that the nuclear region of 
M 33 is on average nearly as old as the outer disc, which they 
found to have an age of ~ 6 Gyr. The peak at 6 Gyr may 
not be real, though: the Mass Pulsation relation suggests 
that low-mass stars do not pulsate much. While this is only 
true in a relative sense (as they live long), the exact level 
and duration of pulsation for those stars may be somewhat 
uncertain. The SFR may thus well have been as high, or 
higher earlier than 6 Gyr ago. We are confident, though, that 
the SFR has dropped significantly since 5 Gyr ago (logt < 
9.7). 

A more recent epoch of enhanced star formation is seen 
too, occurring ~ 200-300 Myr ago (logt = 8.3-8.5), nearly 
reaching the same level as the 6-Gyr peak (log t = 9.8). Since 
then the rate of star formation appears to have diminished 
again. This would confirm the early recognition by Minniti 
et al. (1993) of a substantial star formation episode in the 
"bulge" of M33 within the past Gyr (logt < 9). It would 
be difficult to reliably recover such short-duration enhance- 
ments in SFR at intermediate ages more than a Gyr ago. 
For instance the higher bin around 1.5 Gyr ago (logt = 9.2) 
might well correspond to such period of enhanced star for- 
mation, but on the basis of the star number statistics (repre- 
sented by the vertical errorbars in Fig. 7) we must conclude 
that such assertion would be unfounded. 

The present-day SFR in the inner disc of M 33 (over a 4 
kpc^ area) has been estimated by Wilson et al. (1991) on the 
basis of Ha and far-IR emission at ^ = 0.01 Mq yr~^ kpc~^. 
At ^ = 0.0056 ± 0.0008 Mq yr"^ kpc"^ both the ancient 
and recent epochs of enhanced star formation that we have 



identified appear to be similar in intensity to the estimates 
by Wilson et al. Theirs are valid for the most recent few 
10^ yr, however. In our sample of variable stars we count six 
stars with ages t < 20 Myr {\ogt < 7.3, birth masses > 12 
M0). The SFR derived from this recently- formed population 
is ^ = 0.007 ± 0.003 Mq yr"^ kpc"^: while the errorbar on 
this value is rather large the value is in excellent agreement 
with the value derived by Wilson et al. on the basis of a very 
different method. This is very reassuring, and to a certain 
degree validates our methodology. 

Based on our SFH, the total mass in stars formed over 
the history of M 33 amounts to 3.9 x lO'^ M© kpc~^, of which 
3.1 X 10^ Mq kpc"^ (or 80%) was formed at ages t > 4 
Gyr (logt > 9.6). In contrast, the most recent epoch of 
enhanced star formation, at t < 500 Myr (logt < 8.7), has 
only produced 2.4 x 10® Mq kpc~^ (or 6%). These estimates 
assume that there has not been a net migration into or out 
of the central square kpc region that we have monitored. 
Especially for the older stars this may not be true, and they 
could have preferentially migrated out of the central region 
as a result of dynamical relaxation (cf. van Loon et al. 2003) . 
So the star formation rate - and total stellar mass formed - 
that we derive for ages t^5 Gyr (logt^9.7) may have been 
underestimated. 

Long, Charles & Dubus (2002) determined that the nu- 
cleus (central few arcsec) of M33 underwent two discrete 
epochs of star formation, one a Gyr ago (logt — 9) which 
formed 4% of the 2 x lO*' Mq in stellar mass in the nuclear 
star cluster (cf. Kormendy & McClure 1993), and another 
only ~ 40 Myr ago (log t ~ 7.6) which added just half a per 
cent to the cluster's mass. Curiously, while timed differently 
(though with a degree of uncertainty), this resembles the in- 
crement in stellar mass (6%) that wo found to have formed 
within the past 0.5 Gyr (logi < 8.7) in the surrounding cen- 
tral disc. At the star formation rate that we derived for the 
past 20 Myr (logf < 7.3), the most recent of the two star 
formation events in the rmcleus would have lasted approxi- 
mately one Myr (log At — 6), which is a typical timescalc for 
star cluster formation. These comparisons are suggestive of 
some connection between the nucleus and surrounding disc, 
and possibly of a more similar timing of events. 



4.1 Variation in assumptions 

In this section we examine how the derived SFH changes if 
we change our assumptions: the selection or omission of red- 
dened stars, the amplitude threshold, the distance modulus, 
and the treatment of the super- AGB stars. 



4-1.1 Reddened stars 

We can examine the contribution of the reddened stars to 
the SFH that wo derive from the large-amplitude variables. 
Our catalogue of variable stars in the central square kpc of 
M33 contains 318 stars that have J — K > 1.5 mag (but 
one of these is removed because it was affected by blend- 
ing), which corresponds to the significant fraction of 39% 
of all variable stars. Those stars need a correction to their 
K-band magnitudes. Much redder stars, though rarer, may 
have been missed entirely by less IR-sensitive surveys. Of 
all reddened stars, 45% were identified as variable. Of the 
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Figure 8. Contribution to the derived star formation history by 
reddened variable stars with J — K > 1.5 mag (in red), compared 
to the original result (in black). 



remaining 55%, some might well also be variable stars not 
identified as such in our survey, though it is likely that some 
fraction of this concerns non-variable stars seen behind in- 
terstellar dust clouds. 

Fig. 8 shows the SFH we derive based solely on reddened 
{J — K > 1.5 mag) variable stars identified in our survey. As 
before, a mass-dependent de-reddening correction has been 
applied to account for a mixture of carbon and oxygen-rich 
stars. Two things stand out immediately: the major his- 
torical star formation epoch at 4-10 Gyr (logt = 9.6-10) 
disappeared and also the SFR in the most recent 250 Myr 
(log t < 8.4) has dropped. On the other hand, reddened stars 
dominate the variable stars associated with the SFR derived 
at intermediate ages, 0.25-2.5 Gyr ago (logt < 8.4-9.4). 
Most of these are expected to be carbon stars. Clearly, mas- 
sive and low-mass stars do not feature prominently among 
the reddened stars. We argue, as we did in Paper I, that the 
low-mass stars may become very red (and thus be missed 
even by our IR survey) but possibly only very late in their 
pulsating phase, whilst massive stars do not envelope them- 
selves in optically-thick dust envelopes except in the most 
extreme, and thus rare cases (cf. van Loon et al. 1997). In 
conclusion, to obtain a complete SFH based on variable stars 
one must account for reddened stars so as to not underesti- 
mate the SFR at intermediate ages, while the earliest and 
most recent SFR is determined by optically bright (er) vari- 
able stars. 

Our decision to only correct for reddening stars with 
J — K > 1.5 mag, and to correct them to a locus of 
J — K — 1.25 mag, may introduce some uncertainty. If ei- 
ther choice were made differently, the corrected brightness 
and hence inferred birth mass and thus age would be difler- 
ent. We had introduced the 0.25 mag "buffer" between the 
locus of unreddened stars and the threshold for correcting 
for reddening, to allow for photometric uncertainty which, 
if corrected as if it were due to reddening would lead to 



systematic deviations. The photometric uncertainty in the 
colour is typically ~ 0.1 mag (Fig. 1). The amplitudes of 
variability are typically 0.6 mag; while the colours vary less, 
with only a few measurements in the bands other than the 
K-band the time-average colours may be inaccurate by a 
few tenths of a magnitude. This, and the morphology in the 
colour-magnitude diagram (Fig. 1), supports a. J — K > 1.5 
mag threshold. One could argue for a shift of this, or of the 
unreddened locus, by 0.1 mag or so, but this would only lead 
to a shift (largely systematic) in unreddened K-band bright- 
ness by a similar magnitude. This would resemble varying 
the distance modulus by such amount. In Section 4.1.3 we 
investigate such shift but by as much as 0.4 mag. The effect 
of a 0.1 mag adjustment appears to be of little significance. 



4.1.2 Amplitude threshold 

How sensitive is the derived SFH to the exact threshold of 
our survey to the detection of variability? To gain some idea 
of this effect, we set a selection threshold on the estimated K- 
band amplitude (cf. Paper I) of > 0.5 mag. This leaves us 
with 581 variable stars (excluding 9 stars affected by blend- 
ing - see Paper I), still nearly 3/4 of the original selection. It 
is likely that more stars were removed that have either low 
or high mass, as the former may not pulsate as vigorously 
(McDonald et al. 2010) while the latter do but not in terms 
of magnitude - which is a relative quantity, in relation to 
their high luminosity (cf. van Loon et al. 2008). 

Fig. 9 shows the resulting SFH. The overall trend is lit- 
tle affected. Apart from an overall modest drop in the SFR, 
the effect seems to affect mostly the massive stars formed 
in the past 400 Myr (logt < 8.6). The smaller amplitudes 
of low-mass stars is not apparent, but our SFH also does 
not extend into the past as far as the ages of old Galactic 
globular clusters (i > 10 Gyr i.e. logt > 10) upon which the 
evidence for diminished pulsation in low-mass red giants is 
largely based. 



4.1.3 Distance modulus 

In our analysis we adopted a distance modulus of = 24.9 
mag, which is found by several authors via different meth- 
ods including eclipsing binaries, and which was found to 
match the theoretical isochrones quite well to our IR colour- 
magnitude diagrams. However, some authors have found 
lower values for the distance modulus, around n = 24.5 mag 
(Scowcroft et al. 2009). Changing the distance modulus will 
change the absolute magnitude of stars and hence also the 
mass of stars, but the Mass-Pulsation relation will remain 
unchanged as it does not depend on distance. This means 
that changing the distance modulus will affect the SFH we 
derive. 

Fig. 10 shows the SFH for fi — 24.9 (as before) and /i = 
24.5 mag. Overall the SFH has a similar behaviour despite 
the 20% difference in distance, with a major ancient epoch of 
star formation and enhanced star formation in recent times. 
However, not surprisingly, the SFH is shifted towards earlier 
epochs. The SFR seems to be diminished at all later times; 
this is caused by the same effect, stars having been shifted 
to earlier times, piling up at ages t 10 Gyr (logt ~ 10). 
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Figure 9. Star formation history derived after imposing a higher 
threshold on the estimated ampHtude of variabiUty, Ak > 0.5 
mag (in red), compared to the original result (in black). 
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Figure 11. Star formation history compared for four choices of 
metallicity; Z = 0.004 (blue), Z = 0.008 (green), Z = 0.015 
(black, appropriate for the M33 core) and Z = 0.019 (red). 
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Figure 10. Star formation history compared for two choices of 
the distance modulus: fj, = 24.9 (black) and fi = 24.5 mag (red). 



4-1-4 Metallicity 

While the metallicity of the pulsating red giants in the cen- 
tre of M 33 is quite reliably known, the oldest pulsating 
stars might be somewhat deficient in metals compared to 
the young and intermediate-age Z = 0.015 population that 
formed after the initial period of star formation from which 
most chemical enrichment resulted. Investigating the metal- 
licity dependence of our results also lays the foundation for 



the application of the method to the outskirts of the M 33 
disc, or to altogether different galaxies. 

The parameterisations of the Mass-Luminosity relation 
(Figs. 2 & 3; Table 1), Mass-Pulsation relation (Fig. 5, top; 
Table 2) and Mass-Age relation (Fig. 6; Table 3) have all 
been determined for four different values of the metallicity: 
Z — 0.004 (appropriate for the Small Magellanic Cloud), 
Z — 0.008 (appropriate for the Large Magellanic Cloud), 
Z — 0.015 (appropriate for the centre of M 33) and Z = 
0.019 (super-solar). The resulting star formation histories 
are shown in Fig. 11. 

The result from adopting Z — 0.008 hardly differs from 
the one we obtained for Z — 0.015. Apart from pushing 
the early epoch of star formation back by one or two Gyr 
there are no important deviations outside a ~ 20% mar- 
gin. Adopting Z = 0.019 mainly extends the early epoch of 
star formation to as recent as 2 Gyr ago (logt — 9.3) - the 
star formation rate would become fairly constant over the 
entire history of M 33 except for a lull in activity ~ 0.4-2 
Gyr ago {logt — 8.6-9.3). This strange result can be dis- 
missed as such high metallicity is an unrealistic choice for 
M 33. Adopting Z = 0.004 yields a "burst" of star formation 
0.2-1 Gyr ago (logt = 8.3-9); whilst not impossible this is 
attributable to the shorter duration of pulsation predicted 
by the models (Fig. 5) for the birth mass range of 2-4 M0 
(logM = 0.3-0.6; these are carbon stars) in combination 
with the observed high numbers of pulsating stars in that 
mass range. This low metallicity is not characteristic of the 
M 33 population. In conclusion, the results we have derived 
by adopting Z — 0.015 are robust against reasonable devia- 
tions from this value for the metallicity. 
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1.5 




Figure 12. Mass— Luminosity relation as in Fig. 3 but with 
the revised parametrisation for the super-AGB stars and low- 
mass red supergiants. The solid lines are linear spline fits, re- 
quiring monotony everywhere except for the discontinuity at 
log(M/MQ) = 0.7. 



4.1.5 Super-AGB stars 

Previously, we have interpolated the Mass-Luminosity and 
Mass-Pulsation relations across the super-AGB and low- 
mass red supergiant regime, as the theoretical models had 
not been computed until the very end of these stars' evo- 
lution and we may thus have missed the coolest, strongly- 
pulsating phases that these stars experience before ending 
their lives. Here, we examine how our results change if we 
do take the endpoints of the theoretical isochrones at face 
value. 

Fig. 12 shows the revised Mass-Luminosity relation: the 
AGB portion, at log(M/M0) < 0.7, remains unchanged but 
the description of the more massive portion has been altered 
to follow the isochrones more closely (Table 4). We still keep 
that portion in itself monotonic, for simplicity. The small 
ranges where the slope has the opposite sign will not make a 
very different contribution because of this. Clearly, however, 
stars oi K 15^^ magnitude have two solutions for the mass, 
and we will need to devise some procedure to distribute the 
stars over these two solutions. 

In the overlap region, we let us be guided by the IMF 
as to the relative fractions of the star that are to be assigned 
to the low- and high-mass solutions, respectively: 



Plow 



Phigh 



/iMF dm 



/imp dm 



(11) 



(12) 



Here, mi and m2 are the masses calculated using the low- 
mass solution for each star, between K — 5K and K + 5K, 
where SK is the error on the K-band magnitude; likewise, 
m'l and m'2 delimit the high-mass range. Hence we obtain the 



Table 4. Mass-Luminosity relation, logM [Mg] = aK + b, where 
K is the K-band magnitude on the 2MASS system, for four values 
of metallicity, for the case in which the endpoints of the isochrones 
covering the super-AGB phase have been taken at face value. 
Only the high-mass portion is given here - the low-mass portion 
remains as in Table 1. 



Z = 0.019 


a 


b 


validity range 


-1.749 


25.034 


K < 13.590 


-0.173 


3.609 


13.590 < K < 15.762 


-0.889 


14.893 


15.762 < K < 15.948 


Z = 0.015 


a 


b 


validity range 


-2.198 


31.305 


K < 13.633 


-0.451 


7.491 


13.633 < K < 13.917 


-0.170 


3.582 


13.917 < K < 15.930 


-1.673 


17.311 


15.930 < K < 16.023 


Z = 0.008 


a 


b 


validity range 


-0.888 


13.486 


K < 13.732 


-0.223 ±0.004 


4.355 ±0.008 


13.732 <K < 14.837 


-0.157 ±0.006 


3.378 ± 0.009 


14.837 < K < 15.942 


-0.623 ±0.003 


10.811 ±0.006 


15.942 <K < 16.211 


Z = 0.004 


a 


b 


validity range 


-0.272 ±0.012 


5.148 ±0.018 


K < 14.453 


-0.681 ±0.012 


11.005 ±0.018 


14.453 < K < 15.288 


-0.183 ±0.013 


3.801 ±0.019 


15.288 <K < 16.124 


-0.458 ± 0.007 


8.222 ±0.014 


16.124 < K < 16.386 



normalised probabilities for the low- and high-mass solutions 
as follows: 



-flow — 



Plo 



-Phigh 



Plow + Phigh 

_ Phigh 

Plow + Phigh 



(13) 



(14) 



so that each star still has unit weight, Piow + fhigh = 1- 

The remaining procedure follows that described above, 
but using the modified Mass-Pulsation relation given in Ta- 
ble 5 and displayed in Fig. 5. The resulting SFH is shown in 
Fig. 13. 

In the overlap region we notice a sudden increase in the 
SFR, to ~ 20 times larger than the levels derived by inter- 
polation. Such extreme behaviour right where the modifica- 
tion is made is highly implausible. Could it result from the 
implicit assumption that the SFR was the same for both so- 
lutions for the star's mass, underlying the formulae for the 
fractions piow and Phigh? The inflated values for the SFR are 
due to the high-mass branch of the dual solution. If the SFR 
in reality declined with time, then smaller fractions of stars 
ought to be assigned to that solution and the SFR would 
come down. However, the SFR must decline precipitously in 
order for the high-mass branch to be depleted, essentially re- 
sulting in fewer than one star to be observed in that regime. 
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Table 5. Mass— Pulsation relation, log((5t/t) = D + 
a^exp ((log Af [Mg] -bi)^/{2cf)y where St is the pulsa- 
tion duration and t the age of the star, for four values of metallic- 
ity, for the case in which the endpoints of the isochrones covering 
the super-AGB phase have been taken at face value. 



Z = 0.019 


D 


i 


li 






validity range 


-4.63 


1 


3.45 


1.273 


0.262 


log M > 0.68 




2 


0.27 


0.611 


0.016 


log M > 0.68 


-4.45 


1 


2.10 


0.538 


0.166 


logM < 0.68 




2 


0.75 


0.149 


0.153 


log M < 0.68 


Z = 0.015 


D 


i 




b. 




validity range 


-4.52 


1 


1.28 


1.840 


0.362 


log M > 0.68 




2 


2.88 


1.211 


0.215 


log AI > 0.68 


-3.72 


1 


1.63 


0.468 


0.180 


log M < 0.68 




2 


-0.53 


0.361 


0.085 


logM < 0.68 


Z = 0.008 


D 


i 


ai 


h. 


Ci 


validity range 


-7.13 


1 


5.66 


1.306 


0.381 


log M > 0.68 


-4.02 


1 


1.62 


0.509 


0.217 


log M < 0.68 




2 


0.30 


0.137 


0.054 


log M < 0.68 


Z = 0.004 


D 


i 


ai 




Ci 


validity range 


-4.11 


1 


2.52 


1.313 


0.195 


log M > 0.68 


-4.44 


1 


1.82 


0.602 


0.142 


log M < 0.68 




2 


1.33 


0.202 


0.132 


log M < 0.68 



- SFR, inLerpolaLion 

= SFR, two-mass solution 



_ii=flnzI_L. 



log t (yr) 

Figure 13. Star formation history for Z = 0.015, derived by 
assigning stars to dual solutions for their mass (in red), compared 
to the original result (in black); the vertical dotted lines indicate 
the extent of the overlap region. 



This would be very coincidental and thus also highly un- 
likely. Rather, we find it more likely that the high values for 
the SFR result from unrealistically-brief adopted durations 
for their variability, St (cf. Fig. 5). 

We note that the older of the three age bins in the 
overlap region (at logt ~ 8.4) now has a somewhat reduced 
SFR: this is because these stars have been partly assigned 
to the high-mass branch of the solution whereas in the inter- 
polated case they were all assigned to the low-mass branch 
of the solution. This is not odd in itself. 

Thus we conclude that consistent use of the isochrones 
throughout the super-AGB and low-mass red supergiants 
regime leads to erroneous results. Thereby validating our 
interpolation of the Mass-Luminosity and Mass-Pulsation 
relations assuming that those stars do become cooler and 
more luminous (cf. Ritossa, Garci'a-Berro & Iben 1996; Siess 
2010) and that they do pulsate over longer durations than 
suggested by the isochrones (cf. Ventura & D'Antona 2011). 
As a consequence, we expect super-AGB stars to develop 
strong, dusty winds (following the mass-loss prescription of 
van Loon et al. 2005). This lends support to the suggestion 
by Botticella et al. (2009) that the dust-enshrouded progen- 
itor of SN 2008S was a super-AGB star. 

We remind the reader that while the super-AGB stars 
undergo HBB and thus reach higher luminosities than cal- 
culated in the isochrones, red supergiants do not experience 
HBB and it is therefore not obvious that their evolution was 
also foreshortened in the calculation of the isochrones. The 
exact mass range of super-AGB stars is uncertain and could 
extend up to as high as ~ 11 M0 if convective overshoot is 
ineffective (Eldridge & Tout 2004). 

4.2 Spatial variations 

Up till now, we have considered the SFH derived from the en- 
tire central square kpc (4' x 4') of M 33. However, this region 
may well show different galactic structures associated with 
populations of stars formed at different times. Figure 1 in 
Paper I clearly shows a central concentration of stars giving 
the appearance of a "bulge" of ~ 1' diameter. To examine 
this, first we define a circular region of radius R — 110", to 
avoid problems caused by the corners of the rectangular area 
covered in our UIST monitoring survey. Secondly, we sepa- 
rate the stars in our catalogue into massive stars, AGB stars, 
and RGB stars, on the basis of K-band magnitude and J-K 
colour criteria (see Fig. 1): we define a demarcation line be- 
tween hot massive stars and cooler, less-massive giant stars 
to run from {J-K, K) = (0.6, 18) to {J-K, K) = (0.9, 15.6) 
mag, with massive stars those that have J-K colours bluer 
than this (down to K = 19.5 mag) or that have K < 15.6 
mag, and AGB stars and RGB stars those that have J-K 
colours redder than this and that have 16 < K < 18 mag 
(AGB) or 18.3 < K < 19.5 mag (RGB, which does extend 
to much fainter magnitudes than our survey completeness). 
We left small gaps to avoid contamination. 

We consider two types of galactic structures: disc- 
related and spheroidal. The former have axi-symmetry in 
the galaxy plane and are otherwise "fiat" (with relatively 
little depth), while the latter have a surface density that is 
invariant to inclination with respect to the plane of the sky 
(as long as internal extinction is negligible, which we have 
shown in Paper I is largely the case for our near-IR survey) . 
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Figure 14. Radial distributions of massive stars (in red, mul- 
tiplied by 3 for ease of comparison), AGB stars (in black), and 
RGB stars (in blue), in the galaxy plane. Each of these are fitted 
with a Sersic profile (the fit parameters are listed in Table 6). 
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Figure 15. Radial distributions of massive stars (in red, mul- 
tiplied by 3 for ease of comparison), AGB stars (in black), and 
RGB stars (in blue), in the image plane. Each of these are fitted 
with a King profile (the fit parameters are listed in Table 6). 



To investigate the disc-related structures, we must depro- 
ject the positions of stars on the image to the galaxy plane 
- but for the spheroidal structures we must make sure not 
to do that. To deproject, we first rotate the image according 
to the position angle of PA — 23° and then stretch the co- 
ordinates along the direction of inclination according to an 
incUnation angle of i = 56° (Zaritsky et al. 1989). 

The radial distributions of stars as measured in the 
galaxy plane are shown in Fig. 14. These are fitted with 
an adaptation for the stellar number density of the surface 
brightness function from Sersic (1963): 



N{r) = No exp (-6„(r/i?o)'''" - l) , 



(15) 



where n = 4 returns the de Vaucouleurs' profile typical of 
galaxy bulges (de Vaucouleurs 1953) and n = 1 an exponen- 
tial disc. The calibration of the Sersic profile is such that 
half the intensity (number of stars in our case) falls within 
r = Re- The parameter b„ can be derived exactly from com- 
parison of the complete and incomplete Gamma functions 
(Ciotti 1991), and approximately from b„ = 1.9992n-0.3271 
for 0.5 < n < 10 (Capaccioli 1989). The fit parameters are 
listed in Table 6. 

The AGB distribution shows clear signs of a double- 
component profile, with the break occurring around r ~ 0.4 
kpc, so in this case we fitted a Sersic profile, with — 0.30 
kpc, plus a pure exponential disc profile: 



N{r) = iVoexp(-r/7?e) . 



(16) 



This corroborates the result found by Kormendy & Kenni- 
cutt (2004) for the surface brightness profile: a Sersic profile 
with n = 1.09 ± 0.18 and Re = 0.31 ± 0.05 dominating the 
inner part with the outer part described by an exponential 
disc. Clearly, the AGB stars do not constitute a "bulge", 
which would have had n closer to 4. Apart from the pure 



(rotating) disc component the nuclear part could possibly 
resemble a pseudo-bulge (which is a disc phenomenon and 
as such has n also close to 1). 

The RGB population shows no strong central conden- 
sation, and a large effective radius i?c ~ 2 kpc. This is easy 
to understand as the RGB stars have had many Gyr to re- 
lax within the overall gravitational potential of the entire 
M 33 galaxy. The massive star population too, while more 
centrally concentrated than the RGB population, is spread 
throughout the M 33 disc, without a strong nuclear compo- 
nent. This suggests that steady star formation takes place 
throughout the M 33 disc, with only a marginal concentra- 
tion towards the nucleus. Possibly, the epochs of enhanced 
star formation are a feature of the central regions of M 33, 
leading to the lasting imprint in the AGB stars' distribution 
in the form of a pseudo-bulge. 

The radial distributions of stars as measured in the im- 
age plane are shown in Fig. 15. These are fitted with one of 
the various King models (King 1962): 



N{r) = iVo (1 + {r/R,y 



(17) 



The fit parameters are listed in Table 6. 

Interestingly, the fits are not significantly better or 
worse than for the Sersic profiles fitted in the galaxy plane. 
While it may seem reasonable to expect that the old RGB 
stars constitute a pressure-supported, more spheroidal sys- 
tem, the success of the Sersic fit with a small value for n 
suggests that perhaps they form part of a more flattened 
system. In either case the RGB population has the largest 
characteristic radius among the three populations. The mas- 
sive stars and at least a fraction of the intermediate-age AGB 
stars are not expected to reside in a spheroidal configuration, 
as they form within the gas-rich disc and will not have had 
time to reconfigure into a dynamically relaxed system - for 
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Table 6. Fit parameters of Sersic and King profiles for mas- 
sive stars, AGB stars (Sersic+disc), and RGB stars in the central 
square kpc of M 33. 





Sersic profile 








Population 


No 








n 




(arcsec^^) 


(pc-2) 


(") 


(kpc) 




Massive stars 


0.0021 


0.00010 


199 


0.92 


0.67 


AGB stars 


0.025 


0.0012 


64 


0.30 


0.78 


(exp. disc:) 


0.050 


0.0023 


846 


3.92 




RGB stars 


0.042 


0.0020 


420 


1.94 


0.52 




King profile 








Population 


No 






Rc 


/3 




(arcsec~^) 


(pc-2) 


(") 


(kpc) 




Massive stars 


0.010 


0.0005 


172 


0.80 


1 


AGB stars 


0.24 


0.011 


47 


0.22 


0.4 


RGB stars 


0.17 


0.008 


229 


1.06 


1 
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Figure 16. Ratio of stellar density inside and outside of 120 pc 
(top panel), 250 pc (middle panel), or 500 pc (bottom panel) for 
the massive stars, AGB stars, and RGB stars. While the RGB 
populations are not centrally concentrated (a ratio near unity) 
the AGB and massive stellar populations are progressively more 
centrally concentrated the younger they are. 



Figure 17. Azimuthal distributions of the massive stars (in red), 
AGB stars (in black), and RGB stars (in blue), in the image plane 
(top panel) and galaxy plane (bottom panel). 



One could wonder whether we are aheady noticing the 
gravitational potential responsible for the spiral arm pattern 
that so clearly defines M 33. Or, would the pseudo-bulge be 
tri-axial, like a bar potential? In either case we would ex- 
pect an azimuthal modulation with two max;ima and two 
minima ofTset by an angle Acp — n. In Fig. 17 we plot these 
distributions for the three stellar populations, both in the 
image plane as well as in the galaxy plane. These distribu- 
tions are constructed only for stars within circular areas in 
the image and galaxy planes, respectively, so as not to intro- 
duce spurious modulations due to the corners of the image 
resulting from the generally falling stellar surface density 
with distance to the galaxy centre. 

The AGB population seems not to show any significant 
azimuthal structure in the image plane, but in the galaxy 
plane it does show the modulation expected from a bar or 
two-armed spiral arm pattern, with stellar surface density 
peaks around (^galaxy ~ 100° and 280°, i.e. roughly in the 
East and West directions. The RGB population shows a sim- 
ilar distribution in the galaxy plane, with an enhancement 
also visible in the image plane around (^iimago ~ 260-330°, 
i.e. Westwards. The massive stars, on the other hand, show 
a clear overdensity towards the South, i/iimage ~ 170°. 



these populations the Sersic profile fits yield more plausible 
results. 

Another way to illustrate the central condensation as 
a function of population age is achieved by comparing the 
average stellar surface density within and outside a given ra- 
dius (Fig. 16; the radius is defined in the galaxy plane). This 
suggests that star formation has gradually become more con- 
centrated towards the central regions, or that star formation 
has always been concentrated towards the core with dynam- 
ical relaxation having dispersed stars as time went by. We 
stress that we are looking here at variations within the in- 
ner square kpc, and that further out the situation might well 
differ. 



5 DISCUSSION AND CONCLUSIONS 

The photometric catalogue of Javadi et al. (2011a) was used 
to reconstruct the star formation history and structure of 
the inner square kpc of M33. The numbers and luminosi- 
ties of the pulsating AGB stars and red supergiants were 
converted to star formation rates as a function of look-back 
time, using recent Padova stellar evolution models (Marigo 
et al. 2008). Consistency checks between the predicted lumi- 
nosities and pulsation durations and the derived star forma- 
tion rates lend strong support for these models except that 
the super-AGB stars must also reach high luminosities and 
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Figure 18. (Loft:) Drawing of M33 made by R.J. Mitchell (Rosse 1850). (Middle:) Ultraviolet Imaging Telescope image in blue, against 
an optical image in yellow/red (Landsman et al. 1992). (Right:) UFTI image of the central square kpc of M 33 (the box in the left panels), 
with overlain the massive stars (circles) and the circular area within which we constructed the azimuthal distribution in the image plane. 



pulsate for a considerable amount of time to avoid implau- 
sible values for the inferred star formation rates. This also 
supports the interpretation of dust-enshrouded supernovae 
as being due to super- AGB star progenitors. 

Our conclusions pertaining to the star formation history 
in the centre of M 33 can be summarised as follows: 

• The disc of M33 was built > 6 Gyr ago (logt > 9.8), 
when > 80% of the historic star formation in M 33 occurred. 

• Shorter epochs of enhanced star formation (may) have 
occurred since then; the most recent substantial example 
of such event occurred ~ 250 Myr ago (logf ~ 8.4) and 
contributed < 6% to M 33's historic star formation. 

• A minor star formation event which had been suggested 
to have occurred within the nuclear star cluster within the 
past 40 Myr {\ogt^7.6) can also be discerned in our data 
suggesting this was a more wide-spread phenomenon. 

The formation epoch of M 33 around t > 6 Gyr is in 
line with the peak in the distribution of formation epochs of 
bulgeless galaxies in a standard A-CDM cosmological model, 
around redshift z ~ 1 (Fontanot et al. 2011) corresponding 
to ages t ~ 7 Gyr. We speculate that more recent star forma- 
tion epochs may have been induced by interaction with the 
neighbouring more massive Andromeda spiral galaxy (M 31; 
Richardson et al. 2011). 

Comparison with the historical drawing of M 33 by R.J. 
Mitchell and published in Lord Rosse's records shows a rea- 
sonable correspondence between the footpoints of the spiral- 
arm pattern revealed in the unresolved white hght and the 
density enhancements that we recorded (Fig. 18). The two 
main spiral arms seem to originate around 0imago 90-120° 
and ~ 270° with a Southern spur visible close to the centre. 
The spiral arms may be associated with the overdensities in 
AGB stars and RGB stars, though it would be somewhat co- 
incidental that the signal in the AGB distribution is so per- 
fectly cancelled by the projection on the sky. Also, one might 
not have expected a signal this clear in the distribution of 
RGB stars, if these occupy a more spheroidal, dynamically- 



relaxed distribution. Possibly we are dealing with a bar-like 
feature, which is a disc-related structure and may be con- 
nected to the footpoints of the spiral arms - in fact the Sersic 
index n ^ 0.5-0.8 is smaller than that of a pure exponential 
disc, with galactic bars generally having n ~ 0.5. 

The Southern spur, on the other hand, can be associated 
with the overdensity of massive stars, which are displayed on 
our UFTI image (Fig. 18). This includes a prominent clus- 
ter at {RA,Dec) ~ (l''33™53", +30°38'15"). The Southern 
spur is also prominent in ultraviolet images obtained with 
the Ultraviolet Imaging Telescope on-board space shuttle 
Columbia (Landsman et al. 1992) . Note also that the region 
within 30" is actually rather devoid of massive stars; this 
is also where the pseudo-bulge (or bar) is seen, which there- 
fore seems to be a predominantly intermediate-age feature. 
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